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Abstract 

A nodal analysis method for simulating inertance tube pulse tube refrigerators is introduced. The energy equation, continuity 
equation, momentum equation of gas, energy equation of solid are included in this model. Boundary condition can be easily 
changed to enable the numerical program calculate thermal acoustic engines, inertance tube pulse tube refrigerators, double inlet 
pulse tube refrigerators, and others. Implicit control volume method is used to solve these equations. In order to increase the 
calculation speed, the continuity equation is changed to pressure equation with ideal gas assumption, and merged with momentum 
equation. Then the algebraic equation group from continuity and momentum equation becomes one group. With this numerical 
method, an example calculation of a large scale inertance tube pulse tube refrigerator is shown. 

© 2004 Elsevier Ltd. All rights reserved. 

Keywords: Pulse tube refrigerator; Regenerator; Numerical simulation; Cryocooler; Acoustic engine 


1. Introduction 

An inertance tube pulse tube refrigerator is a very 
important type pulse tube refrigerator [1], There are 
several parameters in this machine. Numerical simula¬ 
tion is an effective way to explain the working mecha¬ 
nism and also an effective way to find the rough 
optimum point for the initial design [2]. In this paper, a 
nodal analysis method based on one-dimensional com¬ 
pressible gas flow model is introduced. The energy 
equation, continuity equation, momentum equation of 
gas and energy equation of solid are included. Acceler¬ 
ation method is used for increasing the calculation 
speed. 

Thermal acoustic engine [3] is a no moving parts 
compressor for pulse tube refrigerator. The development 
of pulse tube refrigerator without moving parts is also 
very important. This method is also designed for cal¬ 
culation of thermal acoustic engine [4] because thermal 
acoustic engine and inertance pulse tube refrigerator are 


*An calculation example is at http://hompage3.nifty.com/cryo- 
cooler/. 
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essentially the same from the numerical viewpoint. It is 
also designed as a general method that can be easily 
changed to simulate other type of refrigerators, such as 
double inlet pulse tube refrigerator, orifice pulse tube 
refrigerator, GM refrigerator, and Stirling refrigerator. 
With this numerical method, a large scale inertance tube 
pulse tube refrigerator is taken as an example calcula¬ 
tion. 


2. Inertance tube pulse tube refrigerator and acoustic 
engine 

Fig. 1 is the schematic of an inertance tube pulse 
tube refrigerator. The volume of compressor 1 changes 
sinusoidal to cause the pressure oscillation within the 
cooler 2, regenerator 3, cold heat exchanger 4 and 
pulse tube 5. The gas in the inertance tube also oscil¬ 
lates. 

Fig. 2 shows an acoustic engine. Heat is supplied to 
the hot heat exchanger 2, and part of the heat is changed 
to work by stack 3, the waste heat is discharged from the 
cold heat exchanger 4. 

From mathematical viewpoint, the above two ma¬ 
chines are the same, only boundary condition is differ¬ 
ent. So the same numerical model is applicable. 
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Nomenclature 


A cross gas flow area 

A s cross solid area 

A l heat transfer area per meter 

c sound velocity 

C v constant volume specific heat of gas 
C P constant pressure specific heat of gas 

C s specific heat of solid 

D diameter of inertance tube 

D e hydrodiameter 

E exergy flow 

/ friction coefficient 

H enthalpy flow 

//px pulse tube enthalpy flow 

// re g regenerator loss 

L length of inertance tube 

L } the distance of the lowest pressure ratio point 
to the left end of the inertance tube 
L 2 the distance of the lowest pressure ratio point 
to the right end of the inertance tube 
LI the last nodal number 

L2 the nodal number in the last control volume 
Lbs the first nodal in the left end of the right buffer 
LPE the last nodal of pulse tube 

LPS the first nodal of pulse tube 

m mass of gas 

m s mass of solid 

m mass flow rate 

P pressure 

Q cooling power 

R gas constant 

5 entropy 

S entropy flow 


1 time 

T temperature of gas 

T 0 room temperature 

7’s temperature of matrix 

u velocity 

V volume 

V 2 volume of compressor 

V 2S dead volume of compressor 

V 20 swept volume of compressor 

V b buffer volume 

V p pulse tube volume 

Vn p i tipi volume 

V tip2 tip2 volume 

W compressor work 

x coordinator along the length 

Greek letters 

6X X step length 

§t time step length 

a heat transfer coefficient 

jl coefficient for equivalent pulse tube volume 

p density of gas 

p s density of matrix 

t period 

a> angle frequency 

Subscripts 

f interface of each control volume 
i number of control volume 

cv center of main control volume 

Superscript 
k time step 



Fig. 1. Inertance tube pulse tube refrigerator. 1. Compressor, 
2. cooler, 3. regenerator, 4. cold heat exchanger, 5. pulse tube, 6. tipi, 
7. inertance tube, 8. tip2 and 9. buffer. 



Fig. 2. Thermal acoustic engine. 1. Hot buffer, 2. heater, 3. stack, 
4. cooler, 5. tipi, 6. inertance tube, 7. tip2 and 8. buffer. 


3. Numerical model 

One-dimensional compressible fluid flow model is 
introduced. The control equations are following [5], 
Continuity equation 


Momentum equation 

Wt {pAu)+ ^ {f>Au2)+A ^c +A D e l 2 PU1 R = ° (1 ' 2) 

Energy equation for gas 


a 

Of 


{p A )+Qx{P uA ) = 0 


0 

( u 2 \ 

0 

/ u 2 \ 

0f 


0X 

puA ( CpT + — J 


+ *A L {T-T t ) = 0 (1-3) 
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Energy equation for matrix 
QT 

AsPsCs-^- = «Al(T - T s ) (1-4) 


4. Numerical method 


Finite difference method is used to solve the above 
equations [6-8]. Each component such as pulse tube, 
regenerator is divided into small control volumes. All of 
the control volumes are put together, then we get the 
mesh system as shown in Fig. 3. Each box means one 
control volume. The black point means the nodal which 
is in the center of the control volume. The temperature, 
mass, density, pressure, and gas property are the values 
at the center of each control volume. The mass flow rate 
and velocity are the values at the interface whose num¬ 
ber is the number of right side nodal number. The vol¬ 
ume between each dotted line is the control volume for 
mass flow rate and velocity. The control volume 2 and 
L2 are changeable in order to change the program easily 
for simulation Stirling refrigerators, GM refrigerators, 
and others. 

In time coordinate, the first order method is used at 
the first step. After step one, the second order method is 
used. In x coordinate, the second order up-wind method 
[9] is used mainly. To avoid the numerical noise, at both 
ends of each component, first order up-wind method is 
used. 

The continuity equation and momentum equation are 
combined. First we change the continuity equation to 
the form P, = /(m,-,m,- + i) under the assumption of ideal 
gas, and combine it with the momentum equation. The 
momentum equation and continuity equation are 
changed to one algebraic equation groups. Compare to 
the method in Ref. [6,7], the algebraic equations groups 
which we have to use the iteration method to solve is 
decreased to one group from two groups. The calcula¬ 
tion speed is increased. 


4.1. Continuity equation 


Eq. (1-1) is changed to the following form with the 
first order method 


k— 1 

rrii — m* 
& 


= m, - m i+ i 


( 2 - 1 ) 


Eq. (2-1) can be changed to the pressure equation with 
the ideal gas state equation by the first order method. 
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(»Z, 


m,-+i)5i + m) 1 



m i+ y)ht + m k i 1 



( 2 - 2 ) 

(2-3) 


Pi = (m, - m,-+i)a,- + b, 


a, = ht 


RTj 

V, 


b, =; 


RT, 

IT 


(2-4) 

(2-5) 

( 2 - 6 ) 


Eq. (1-1) can be changed to the pressure equation by the 
second order method 


1.5m,-2m^ 1 +0.5mf- 2 


5 1 


= m, — m i+ 1 


1.5 V,Pi 
RT: 


= 1.5m, = (/»,• — m i+ i)8f + 2m] — 0.5m' 


Jfc-2 


Pt = 


(m, — m i+ i)bt + 2m * 1 — 0.5 m) 


RTj 

J TJV: 


Pi = (m, - m,+i)fl/ + hi 
s RT, 

a, = 8 t 


1.5 V, 


bi = (2.0m* _1 — 0.5m* -2 ) 


(2-7) 

( 2 - 8 ) 

(2-9) 

( 2 - 10 ) 

( 2 - 11 ) 

( 2 - 12 ) 


4.2. Momentum equation 
Eq. (1-2) can be changed to 

0 , x 0 QP / 1 

+ -mlul = 0 

Qjfl ■ 

-Qf bx i + 0”«)cV/ - 0«)cvi-l +Ai(Pi-Pi- 1) 

8 Xj 1 . . 

+ ~D~f i 2 lit ‘\ Ui \ = ° 


(3-1) 


(3-2) 


(mu) cyi and (mu) CYi _ 1 is on the nodal point. 

Each term in Eq. (3-2) can be changed to the 
following term. 

By the first order method 


i,k -1 


0 m, m, - m" 

By the second order method 


(3-3) 


0m,- 1.5m, — 2m) 1 + 0.5mf 2 

it =-s- 


Fig. 3. Grid and control volume of regenerator. 


(3-4) 
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By Eq. (2-4) or (2-10) 

Ai(Pi - P- 1) = A l a,m l - T,a,fft,- + i + T,-h, - T,<2,-_im,_i 

+ AiOi-itki - (3-5) 

The term (mu) cvi — (»jm) cv ,_| is a little complex. The 
second orders up-wind method has been used except 
both ends of each component. 

/ a u m, - hum,--! m cv , >0 , , 

( a 2i m i+ 1 - b 2i m i+ 2 m cv ,- < 0 


f au-i/Mf-i - bu-irhi- 2 m cv ,_i > 0 
\ « 2 /- 1 rw, - b 2i -\m i+ \ Mcv,.! < 0 


(3-7) 


Mi,- = 1.5, hi,- = 0.5, a 2 ; = 1-5, b 2i = 0.5 for second order 
method, a u = 1, hi, = 0, a 2i = 1, b 2i = 0 for first order 
method. 

The velocity is changed to the following form 


I M cv , M cv / 0 
{ 0 M CVI - < 0 


(3-8) 


Bi,= 


0 Uq v/ ^ 

^cv/ ^cv/ M ^ 


Then 


(3-9) 


(*k)cw ” (" w )cv/-i = Auiau/hi - h h -m,-_i) 

+ 5 i,-(a 2 ,w/+i - b 2 i m i+2 ) 

— Aii-i(a\i-imi-\ - hi,-_im,_ 2 ) 

- fii,_i(a 2 i _im,- - b 2 i -\th i+ \) 

(3-10) 

If the flow is laminar flow with the friction equation 


Re 


(3-11) 


Then 


Sx,- 

D ei 



/hi / f '/tr bx, 

Dip 


(3-12) 


If the flow is turbulent flow, the friction equation can be 
getting from Ref. [10], 

Then the momentum equation can be changed to 

a p ,-/M, = + a ei m i+l + h p , (3-13) 

The momentum equation is not considered in volume 2 
and volume L2. 


By the first order method in time direction 
0 C v ,m, 7 ) _ C y ,m,T, - C k ^m k r l Tf-' 

dt 8t 1 ‘ j 

By the second order method in time direction 

6Cy,w,-T _ 1.5Cy,/ W ,-T - 2.0C(,7 1 mf- 1 ^- 1 + O.SC^-mf 2 ^- 2 

S t 8 1 

(4-3) 

In x direction, the second order up-wind method are 
used except both ends. 

T _ f ayTj-i - b 2i Ti- 2 nii> 0 , . 

f ' \ a 4 ,-7} - h 4 ,7) + i m, <0 


f a 3i+l Ti ~ bn+lTi-i »?,- + i > 0 

\ «4/+i 7/+i — h 4i+ i7} + 2 h?,- + i < 0 


fl 3 ,- = 1.5, h 3; - = 0.5, a 4( - = 1.5, h 4( - = 0.5 for the second 
order method, a 3i = 1, h 3( - = 0, a 4 , = 1, h 4 , = 0 for the 
first order method. To avoid numerical noise, at both 
end of each component such as regenerator, or the in- 
ertance tube, the first order up-wind method and the 
second order up-wind method are mixed. Fig. 4 shows 
how to use the first order up-wind method and the 
second order up-wind method near the two components 
such as the regenerator and the cooler. Fig. 5 shows how 
to use the first order up-wind method and second order 
up-wind method at both end of the machine. In Figs. 4 
and 5, the arrow to right means gas flow to right 
direction, the arrow to left means the gas flows to left 
direction. The 2nd means the second order up-wind, the 
1 st means the first order up-wind method. 

The mass flow rate timing heat capacity is changed to 
the following 


j Cp,m, rhi > 0 
( 0 ihi < 0 


(4-6) 
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Fig. 4. First order up-wind method near the interface of two compo- 
4.3. Energy equation of gas nents. 


Energy equation can be changed to equation 
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Fig. 5. First order up-wind method at both end of the machine. 
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f 0 ihj > 0 

\ C Pl m, ihj < 0 


We can use the same method to deal with the left 
(4-7) buffer. 


Then 

Cpi+\nti+\Tu+i — Cp/W/Tf,- = A 2i+ \ (<23,-+i T t — b-n + iT t -\) 

+ Bii+\(a4i + \T i+ \ — bn i+ \T i+ j) 

— A 2i (ay 7)_i — Z>3/7]_2) 

— Bsi{ci4iTi — b4,7/ + i) 

(4-8) 

The final equation of the energy equations of gas and 
matrix are changed to 

Op, 7; = a ei 7)_i + a wi T i+ i + b pi (4-9) 

4.4. Energy equation of solid 

Matrix equation (1-4) also can be changed to 
-^ + M t ,{7 Sl -f,) = 0 (5-1) 

By the fist order method 

T r* T rk -1 T'k -1 

0C Si-*Sf _ ISMS/ C S|: 7 Si /r 9 X 

0r 5r 1 

By the second order method 
6Cs,-r s ,- _ 1 .5Cs,T s , - 2.0C*n7*-i + 0.5C*- 2 ?*- 2 
0r Sr 1 ’ 

Finally we get 

Op, 7s; = £ p , (5-4) 

4.5. Buffer 


4.6. Boundary 

The volume of the compressor 

V 2 = V 2S + 0.5F 20 (1 - sin (cor)) (7-1) 

V 2 is a part of buffer if it is an acoustic engine. 

The volume of L2 also can be changeable. The pres¬ 
sure drop equation and heat transfer equation are taken 
from reference [10], The wall temperatures of the heat 
exchangers, buffers and inertance tube are constant. 

4.7. Total parameters 


Enthalpy flow 

H = - j) m^Cp + ^w 2 ^dr (8-1) 

Entropy flow 


S = - j> ms At 

(8-2) 

Exergy flow 


E = H-ST 0 

(8-3) 

Compressor input work 


W = - jP 2 AV 2 

(8-4) 

Q = 7/px — 77 reg 

(8-5) 


There is alternate of the dealing with the momentum 
equation in the right buffer. The momentum effect in the 
buffer is almost zero. We can consider the pressure in the 
buffer is the same and the velocity is zero on each time 
step. It is not necessary to solve the momentum equa¬ 
tions in the buffer. The calculation speed will be in¬ 
creased. 

Apply Eq. (2-2) in the buffer from i = Lbs to L2, then 
we get the pressure in the buffer with the first order 
method. 


COP = Q/W (8-6) 

Equivalent PV diagrams are given from the both ends of 
the gas piston within the pulse tube. The gas piston is 
defined as follows. The maximum left position of the left 
end of the gas piston just touch the left end of the pulse 
tube. The maximum right position of the right end of the 
gas piston just touch the right end of the pulse tube. The 
following term is introduced for getting the position of 
the left and right end of the gas position. 


T Lbs — 


(w L bs - »Ul)Sr+Xjff L bs OT f 1 
2^i =Lbs RTi 


( 6 - 1 ) 


Apply Eq. (2-8) in the buffer from i = Lbs to L2, then 
we get the pressure in the buffer with the second order 
method. 


TVhs = 


(m Lbs - m L \)8t + 2£. 


12 


„k -1 


-0.5£L 2 l 


1 5 V 4L 

1 • J Z^;=Lbs RTi 


( 6 - 2 ) 


LPE 

^=<h t +i>? ( § - 7 ) 

i—i 

where m nght is the gas within the imaginary cylinder at 
the hot end of the pulse tube. 

k 

^right = m right0 + ^LPE+l (^'^) 

k= 1 

^righto is chosen to let the minimum of m k righl become 
zero. 
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The left end of the gas piston is represented by 

M lps = min (M pps ) (8-9) 

The right end of the gas piston is represent by 

M lpe = max (M* pe+1 ) (8-10) 

If M\ < M lps < Mf +1 , we can know the left end position 
of the gas piston 

x lps =* + w - m lps ) (8-n) 

lvl i lvl i+\ 

If M- < M\ VV < M* +l , We also can know the right end 
position of the gas piston 

^pe = * + uf-Ml, {M ' ~ MlPe) (8 ‘ 12) 


4.8. The solve process 

The present program cannot let the gas oscillate by 
itself in acoustic engine. An initial small velocity is 
given in the program. The solving processes are as 
follows: 

(1) Temperature field 

(2) Mass flow rate field 

(3) Pressure field 

It is repeated from (1) to (3) until one time layer is 
finished, and then goes to next time layer till one period 
is finished. The calculation is continued till the required 
calculation accuracy that is judged by how near of the 
enthalpy flow rate at both end of regenerator and pulse 
tube. At each iteration step, tri-diagonal matrix algo¬ 
rithm (TDMA) method is used for solving Eqs. (3-13) 
and (4-9). 

For acoustic engine, during the calculation, the 
maximum pressure point in the hot buffer is checked out 
to get the oscillation period, which is given from two 
maximum pressure points. 


5. Other regenerator type refrigerators 

The present simulation program is also applicable for 
other regenerator type refrigerators, such as Stirling type 
pulse tube refrigerators, GM type pulse tube refrigera¬ 
tors, Stirling refrigerators, GM refrigerators, VM 
refrigerators, and so on. The difference is boundary. The 
momentum term is not necessary in this kind of refrig¬ 
erators. 


6. Calculation example of large scale inertance tube pulse 
tube refrigerator 

Here the big pulse tube refrigerator is taken as an 
example whose regenerator with diameter 110 mm and 
length 50 mm and pulse tube with diameter 38 mm and 
length 150 mm is about 20 times large than those in Ref. 
[11], In Ref. [11], the cooling power is less then 20 W, 
but the cooling power is hundred watts in the present 
refrigerator. Thus we call this pulse tube refrigerator as 
a large scale inertance tube pulse tube refrigerator. 
Room temperature is 300 K. Refrigeration temperature 
is 80 K. Working medium is helium gas. Average pres¬ 
sure is 2 MPa. Working frequency is 50 Hz. The wall 
temperature of the compressor, the cooler, the tips, the 
inertance tube, and the buffer is room temperature. The 
wall temperature of the cold heat exchanger is refriger¬ 
ation temperature. Regenerator matrix and pulse tube 
wall only exchanges heat with helium gas. The expan¬ 
sion work of the pulse tube is very large in the big in¬ 
ertance tube pulse tube refrigerator. The gas in the 
inertance tube can strongly oscillate due to the big work 
input to the inertance tube. 

6.1. Optimum point 

By Ref. [2], if other parameters are fixed, there is an 
optimum diameter and length of the inertance tube. Fig. 
6 a-d shows the inertance tube length dependence of the 
input power of the compressor, cooling power, pressure 
ratio in the pulse tube, and COP, for different diameter 
for the inertance tube when the buffer volume is 4 1 and 
the compressor swept volume is 150 cc. The peak values 
of the cooling power, pressure ratio in the pulse tube 
and COP become the largest with the inertance tube 
with 16 mm in diameter. The input power with inertance 
tube of 16 mm in diameter is only slightly lower than 
that with inertance tube of 20 mm in diameter. So the 
inertance tube of 16 mm in diameter is considered as 
optimum. With inertance tube of 16 mm in diameter, the 
COP gets peak of 0.165 with length 3.75 m, and the peak 
cooling power over 600 W and peak input power about 
4000 W gets with length 3.375 m. 

With inertance tube in diameter 16 mm and length 
3.75 m, the equivalent PV diagrams within the pulse 
tube are shown in Fig. 7, the mass flow rate and pressure 
wave at the cold end of the pulse tube and the center of 
the regenerator are shown in Fig. 8. Compared to the 
small pulse tube refrigerator in Ref. [2], a large phase 
lead of the pressure wave is achieved relative to the mass 
flow rate at the cold end of the pulse tube, and the PV 
diagram at the cold end is tilted to the left largely. The 
pressure wave and mass flow rate are almost in phase at 
the middle of the regenerator. So on the optimum point 
of large scale inertance tube pulse tube refrigerator, the 
PV diagram at the cold end of the pulse tube is tilted to 
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Fig. 6. (a) Compressor input work: D12 inertance tube diameter 12 mm, D16 inertance tube diameter 16 mm, D20 inertance tube diameter 12 mm, 
D24 inertance tube diameter 16 mm. (b) Cooling power: D12 inertance tube diameter 12 mm, D16 inertance tube diameter 16 mm, D20 inertance 
tube diameter 12 mm, D24 inertance tube diameter 16 mm. (c) Pressure ratio in pulse tube: D12 inertance tube diameter 12 mm, D16 inertance tube 
diameter 16 mm, D20 inertance tube diameter 12 mm, D24 inertance tube diameter 16 mm. (d) COP: D12 inertance tube diameter 12 mm. D16 
inertance tube diameter 16 mm, D20 inertance tube diameter 12 mm, D24 inertance tube diameter 16 mm. 



Normalized length of pulse tube 

Fig. 7. Equivalent PV diagrams in pulse tube. 


the left, the pressure wave leads a large phase to mass 
flow rate at the cold end of the pulse tube, and mass flow 
rate and pressure wave in phase point is inside of the 
regenerator. 

6.2. Inertance tube 

The inertance tube is an important part in the iner¬ 
tance tube pulse tube refrigerator. The pressure wave, 
mass flow rate and velocity in the inertance tube is useful 
for understanding how the inertance tube working in 
large scale pulse tube refrigerator. 

With inertance tube in diameter 16 mm and length 
3.75 m, the maximum and minimum pressure and 
velocity in the inertance tube are shown in Fig. 9. Flere 
X is normalized by the total length of the inertance tube 
directed from left to right. Also, the mass flow rate and 
pressure at X = 0, 0.5, 0.87, and 1.0 are shown in 
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0 90 180 270 360 


(b) Normalized time, degree 

Fig. 8. (a) Pressure wave and mass flow rate at pulse tube cold end and 
(b) pressure wave and mass flow rate at regenerator center. 

Fig. 10a and b. As shown in Fig. 9, the pressure 
amplitude decreases from the left end to the lowest 
pressure point (X = 0.87) where the pressure amplitude 
is the smallest, and increases again, while the velocity 
shows a maximum near the right. Fig. 10a shows that 
the pressure amplitude decreases from the left to the 
center while the phase does not vary much. At X = 0.87 
frequency is twice as large as the base frequency because 
the amplitude at the second harmonics is larger than 
that at the base frequency. At the right end, the phase of 
the pressure wave is changed by 180° from the others. 
Fig. 10b shows that the phase of the mass flow rate is 
kept almost constant throughout the inertance tube. In 
most part of the inertance tube, the phase angle between 
the mass flow rate and pressure wave is about 90°. So 
when the velocity is at the maximum, the pressure is 
almost at zero amplitude, when the pressure is at the 
maximum amplitude, the velocity is close to zero. Fig. 
11 a and b show the instantaneous pressure and velocity 
distribution in the inertance tube, respectively. Indicated 
degree correspond to the instance of a cycle which is 
divided into 360°. It is shown that the pressure is almost 
linear distribution in the inertance tube in any time. The 


Fig. 9. Maximum and minimum pressure and velocity in inertance 
tube. 




Fig. 10. (a) Pressure wave in inertance tube: 1. hot end of pulse tube, 
2. center of inertance tube. 3. lowest pressure ratio point and 4. buffer, 
(b) Mass flow rate in inertance tube: 1. hot end of pulse tube, 2. center 
of inertance tube, 3. lowest pressure ratio point and 4. buffer. 


peak velocity is about 100 m/s. So large inertance effect 
could be obtained. 
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(a) Normalized length of inertance tube and tips 



(b) Normalized length of inertance tube and tips 

Fig. 11. (a) Pressure distribution in inertance tube, (b) velocity 
distribution in inertance tube. 
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Fig. 12. Exergy flow and enthalpy flow. 



Normalized inertance tube length 

Fig. 13. Equivalent PV diagrams in inertance tube. 


When the inertance tube of 16 mm in diameter and 
3.75 m in length, the exergy flow and enthalpy flow 
along the refrigerator are shown in Fig. 12. In Fig. 12, 
coordinate X for the different components in the 
refrigerator is normalized by their respective length be¬ 
cause the length of the regenerator and pulse tube is 
millimeter order while the length of the inertance tube is 
meter order. The range 0 < X < 1 represents the cooler 
and the regenerator; the range 1 < X < 2 represents the 
cold heat exchanger and the pulse tube; the range 
2 <X < 3 represents the tipi, inertance tube tip2, and 
buffer. The exergy flow near 3000 W is inputted from the 
compressor. It decreases to about 2000 W at the cold 
end of the regenerator. About 1500 W of it is taken out 
as cooling power through the cold head exchanger. 
Finally, the exergy flow over 500 W flows from the left 
end of the pulse tube into the inertance tube and grad¬ 
ually decreases to zero at the end of the inertance tube. 
The over 500 W exergy flow is the expansion work of the 
gas in the cold end of the pulse tube. It is the power for 
the gas oscillation in the inertance tube. Fig. 12 also 
shows that the loss in the regenerator is very large even 
in the large scale inertance tube pulse tube refrigerator. 


Similar to the pulse tube, we also can get equivalent 
PV diagrams in the inertance tube as shown in Fig. 13 
for the gas pistons dividing the inertance tube into three. 
From left end to right end, the size of the PV diagrams 
decreases, the shape is also changed. 

The numerical results also show that the distance 
from the lowest pressure ratio point to the buffer almost 
does not change regardless of the length of the inertance 
tube, but it increases with the increasing of the diameter 
of the inertance tube when the buffer volume is 4 1. In 
other words, if the diameter of the inertance tube and 
the buffer volume is fixed, the distance from the lowest 
pressure ratio point to the hot end of the pulse tube is 
changed if we change the length of inertance tube. 

6.3. Buffer size and compressor size effect 

With the inertance tube of 16 mm in diameter, the 
input power of the compressor, cooling power, pressure 
ratio in the pulse tube, and COP are shown in Fig. 14 for 
different buffer sizes (1, 2, 4, 10 1) and swept volume 
(150, 200, 250 cc) of the compressor. With decreasing 
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(a) 


Inertance tube length, m (b) 


Inertance tube length, m 



(c) Inertance tube length, m 



(d) Inertance tube length, m 


Fig. 14. (a) Buffer volume and compressor swept volume effect to compressor input work: B10C150 buffer volume 10 1 and compressor swept volume 
150 cc, B4C150 buffer volume 4 1 and compressor swept volume 150 cc, B2C150 buffer volume 2 1 and compressor swept volume 150 cc, B1C150 
buffer volume 1 1 and compressor swept volume 150 cc, B4C200 buffer volume 41 and compressor swept volume 200 cc, B4C250 Buffer volume 41 and 
compressor swept volume 250 cc. (b) Buffer volume and compressor swept volume effect to cooling power. B10C150 buffer volume 10 1 and com¬ 
pressor swept volume 150 cc, B4C150 buffer volume 4 1 and compressor swept volume 150 cc, B2C150 buffer volume 2 1 and compressor swept 
volume 150 cc, B1C150 buffer volume 1 1 and compressor swept volume 150 cc, B4C200 Buffer volume 4 1 and compressor swept volume 200 cc, 
B4C250 buffer volume 41 and compressor swept volume 250 cc. (c) Buffer volume and compressor swept volume effect to pressure ratio in pulse tube: 
B10C150 buffer volume 10 1 and compressor swept volume 150 cc, B4C150 buffer volume 4 1 and compressor swept volume 150 cc, B2C150 Buffer 
volume 2 1 and compressor swept volume 150 cc, B1C150 buffer volume 1 1 and compressor swept volume 150 cc, B4C200 Buffer volume 4 1 and 
compressor swept volume 200 cc, B4C250 buffer volume 4 1 and compressor swept volume 250 cc. (d) Buffer volume and compressor swept volume 
effect to COP. B10C150 Buffer volume 10 1 and compressor swept volume 150 cc, B4C150 Buffer volume 4 1 and compressor swept volume 150 cc, 
B2C150 buffer volume 2 1 and compressor swept volume 150 cc. B1C150 buffer volume 1 1 and compressor swept volume 150 cc, B4C200 Buffer 
volume 4 1 and compressor swept volume 200 cc. B4C250 Buffer volume 4 1 and compressor swept volume 250 cc. 


the buffer size, the peak values of the input power and 
cooling power decrease, but the peak COP value re¬ 
mains almost the same. The optimum length of the in¬ 
ertance tube is increased. When the compressor swept 
volume increases, the input power, cooling power and 
the pressure ratio in the pulse tube are increase, while 
the COP decrease, and the optimum length of inertance 
tube is almost not changed. 

With compressor swept volume 150 cc, other 
numerical results show that when the buffer volume is 
10 1, the optimum diameter of inertance tube is 20 mm 
which is just a little better than the inertance tube of 


16 mm in diameter. The optimum diameter of inertance 
tube is 16 mm in diameter with buffer volume of 1 and 

2 1 . 

With buffer volume of 4 1, other numerical results 
show that the compressor swept volume almost has no 
influence on the optimum of inertance tube. 

Fig. 15 shows the axial distribution of the pressure 
ratio in the inertance tube of 16 nun in diameter for 
different buffer sizes (1, 2, 4, 10 1). The inertance tube 
length are 5.0, 4.25, 3.75, and 3.5 m, respectively, which 
are chosen to give the peak COP values for each buffer 
volume. The swept volume of the compressor is 150 cc. 
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Length of inertance tube, m 


Fig. 15. Pressure ratio in inertance tube. BIO buffer volume 10 1, 
B4 buffer volume 4 1, B2 buffer volume 2 1, B1 Buffer volume 1. 

The distance from the lowest pressure ratio point to its 
left end is almost not changed with the decrease of the 
buffer volume, or the lowest pressure point to the right 
end, is increased in spite of the decrease of the buffer 
volume. 

When the buffer volume is changed, the maximum 
velocity point in the inertance tube is also changed. The 
maximum velocity point is near the lowest pressure ratio 
point which is not shown here. 

6.4. Inertance tube and Helmholtz resonator 

The working mechanism can be explained more 
clearly with Helmholtz resonator [12]. In inertance tube 
pulse tube refrigerator, it can be considered there are 
two helmholtz resonators. As shown in Fig. 16, the right 
part of line 0-0 is a helmholtz resonator with the L 2 
length inertance tube as the resonator. The part between 
line 1-1 and 0-0 is another helmholtz resonator with the 
L\ length inertance tube as the resonator and the volume 
of pulse tube and part of the regenerator as the buffer 
with line 1-1 as the imagine bottom. At line 1-1, the 
mass flow rate and pressure wave is in phase. At line 0- 
0 , it can be consider there is an infinite buffer volume. 
The pressure ratio around line 0-0 is the lowest in the 



Fig. 16. Inertance tube and Helmholtz resonator. 1. Compressor, 
2. cooler, 3. regenerator, 4. cold heat exchanger, 5. pulse tube, 6. tipi, 
7. inertance tube, 8. tip2 and 9. buffer. 


inertance tube, while the velocity is the maximum in the 
inertance tube. 

In Fig. 16, if the mass flow rate and pressure wave is 
in phase at line 1-1, the pressure wave must leads a 
phase angle at the cold end of the pulse tube. 

The frequency of the two helmholtz resonator is the 
frequency of the compressor which is fixed by the de¬ 
sign. The length of L\ and L 2 can be estimated from Ref. 
[ 12 ]. 


71 D 2 / c \ 2 

(9-1) 

4 /lip + ktipl 'ffl' 

b 

(9-2) 

4 V b + V tip2 KcoJ 


Here fl l'\> means the equivalent volume of the pulse tube. 
We have chosen f) as 3 based on Ref. [13]. The total 
length of the inertance tube L is equal to L\ plus L 2 . 

Eqs. (9-1) and (9-2) are useful for the explaining of 
the inertance tube working condition. In Fig. 16, if the 
diameter of the inertance tube and L 2 are fixed, the de¬ 
crease of the total length of inertance tube leads to the 
decrease of L u then the imaginary buffer bottom line 1-1 
will move to the left. If the diameter of the inertance 
tube increased, the length of L 2 is increased, the opti¬ 
mum length of L\ also will be increased because the 
imaginary buffer bottom line 1-1 will not be changed so 
much for the optimum COP. 

Eqs. (9-1) and (9-2) are also useful for a rough esti¬ 
mation the inertance tube. This helps us to decrease the 
numerical calculation job because the calculation speed 
is not so high even we use the acceleration method. The 
length of Ei, L 2 , and L determined by the numerical 
simulation (makers) and by Eqs. (9-1) and (9-2) (curves) 
are shown in Figs. 17 and 18. It is shown that Eq. (9-2) 
has a very good accuracy, while Eq. (9-1) gives much 



Inertance tube diameter, mm 

Fig. 17. Inertance tube length and its lowest pressure ratio point vs. 
inertance tube diameter. 
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Buffer volume, litre 

Fig. 18. Inertance tube length and its lowest pressure ratio point vs. 
buffer volume. 

different values from numerical results when the iner¬ 
tance tube diameter is over 16 mm. 

The accuracy of Eq. (9-1) may be depend on the value 
of (0.257tZ) 2 Ti)/(/?fp + Kipi)- With diameter 16 mm, the 
(0.25nD 2 L[)/(pVp + Ki P i) by Eq. (9-1) is less than 1.2. 

The length of L\ should be less than one quarter wave 
length of helium gas of 5 m at 50 Hz and 300 K in order 
to let the line 1-1 in Fig. 16 in the regenerator. So if the 
calculation from Eq. (9-1) is over 5 m, it should be 
chosen less than 5 m. 


7. Conclusion 

A nodal analytical method based on implicit control 
volume concept is introduced for inertance tube pulse 
tube refrigerator. It can simulate the inertance tube pulse 
tube refrigerator and the thermal acoustic engine. It be¬ 
comes a general numerical method, which can be used 
for the simulation of the double inlet pulse tube refrig¬ 
erator, GM refrigerator, Stirling refrigerator, and other 
type refrigerators, too. The continuity equation is 
changed to pressure equation and merged with momen¬ 


tum equation for decreasing the algebraic equations to 
one group for increasing the calculation speed. A large 
scale inertance tube pulse tube refrigerator is taken as an 
example for the calculation. 
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